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Tidal disruption rate of stars by supermassive black holes 
obtained by direct N-body simulations 
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ABSTRACT 

The disruption rate of stars by supermassive black holes (SMBHs) is calculated nu- 
merically with a modified version of Aarseth's NB0DY6 code. Equal-mass systems 
without primordial binaries are treated. The initial stellar distribution around the 
SMBH follows a Sersic n = 4 profile representing bulges of late type galaxies as well 
of early type galaxies without central light deficits, i.e. without cores. In order to infer 
relaxation driven effects and to increase the statistical significance, a very large set of 
N-body integrations with different particle numbers N, ranging from 10 3 to 0.5 ■ 10 6 
particles, is performed. Three different black hole capture radii are taken into account, 
enabling us to scale these results to a broad range of astrophysical systems with relax- 
ation times shorter than one Hubble time, i.e. for SMBHs up to M. « 10 7 M Q . The 
computed number of disrupted stars are driven by diffusion in angular momentum 
space into the loss cone of the black hole and the rate scales with the total number 
of particles as 4j oc N b : where b is as large as 0.83. This is significantly steeper than 

the expected scaling 4j oc In (AT) derived from simplest energy relaxation arguments. 
Only a relatively modest dependence of the tidal disruption rate on the mass of the 
SMBH is found and we discuss our results in the context of the M, — a relation. The 
number of disrupted stars contribute a significant part to the mass growth of black 
holes in the lower mass range as long as a significant part of the stellar mass becomes 
swallowed by the SMBH. This also bears direct consequences for the search and exis- 
tence of IMBHs in globular clusters. For SMBHs similar to the galactic center black 
hole Sgr A*, a tidal disruption rate of 55 ± 27 events per Myr is deduced. Finally 
relaxation driven stellar feeding can not account for the masses of massive black holes 
M. > 10 7 M© in complete agreement with conventional gas accretion and feedback 
models. 

Key words: black hole physics, spherical galaxies, Sersic profiles, methods: N-body 
simulations, gravitational dynamics 



1 INTRODUCTION 

The evolution of supermassive black holes (SMBHs) 
and their host galaxies is at present one of the key 
problems of astrophysics. Motivated by empirically 
found scaling relations between properties of galaxies 
in terms of velocity dispersion a |Gebhardt et al.l l200d ; 
iFerrarese fc Merrittll2000l ; iGiiltekin et alJ^OOSn T luminosity 
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L l|Kormendv fc Richstonel Il995l; IFerrarese fc Fordl 120051). 



bulge mass MB U i gc ( Magorrian et al. 1998; lHaring fc Rbd 
2004). central ligh t defi ci t L def (lLauer et al l 120071 : 



Kormendv fc Bended 120091 ; iHopkins fc HernquistJ |201Ch 



total number of globular clusters Agc ( Burkert fc Tremaind 
l2010t ) and the mass of their central black holes M. , there is 
a substantial need to understand the related evolution of 
both SMBHs and their hosts. In order to constrain galaxy 
formation models and to answer the question as to what 
powers the growth of SMBHs over cosmic times, all forms 
of matter which are accreted must be taken into account. 
This becomes more urgent as recent studies have found 
evidence for deviations from the general scaling relations 
for the most-massive and for the least-massive black holes 
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dLauer et alj 120071 ; iGebhardt et al.1 l201ll ; iKormendv et all 



Gas accretion is thou ght to be th e most dominant 
driver of SMBH gro wth ijsoltanl Il982ft . Modern studies 
|Yu fe Tremaindl2002T l estimate the black hole mass density 
from the spatial distribution and from the measured stellar 
velocity dispersions in elliptical galaxies in combination 
with the M, — a relation. The SMBH mass density is then 
compared with the observed quasar luminosity function 
in order to yield constraints on the accretion efficiency 
parameter e as well as on the growth history. In order 
to make these studies even more accurate, the impact 
of other feeding modes like merging supermassive black 
holes and stellar captures must also be taken into ac- 
count. Simultaneously the luminous gas accretion history 
of low-mass SMBHs (Af. « 10 s - 10 7 M Q ) is harder to 
measure especially at large redshifts as they never approach 
luminosities comparable to those of quasars. It is even 
plausible that low-mas s SMBHs gain most of th eir mass by 
tidal disruption events |Milosavlievic et al.ll2006h . Therefore 
it is important to infer the stellar capture rate for as many 
astrophysical systems of interest as possible, for all relevant 
SMBH masses using both theoretical and when possible 
numerical approaches. In order to avoid confusion regarding 
the terminology of the capture and disruption rate we note 
that the former expression is used for the general number 
of stars/particles which are either swallowed as a whole 
or disrupted outside the event horizon in a given time, 
i.e. independent of the mass of the SMBH. The latter one 
is explicitly used for situations in which stars are tidally 
disrupted before they would enter the event horizon. 

In this paper we present the disruption rate of stars by 
SMBHs with masses in the lower range up to Mm < 10 7 Mq 
embedded inside realistic stellar density profiles. These 
results are obtained by self-consistent direct N-body 
integrations and increase the hitherto probed region of 
direct numerically inferre d disrup tion rates. Pioneered by 
iBaumgardt et all (|2004al lbl, |2006h for intermediate-mass 
black holes (IMBHs) at the centers of globular clusters, 
our calculations can be applied to a larger sample of 
systems. Our findings should be regarded as comple- 
mentary to other cont r ibutio n s llDuncan fc Shapiro! 1 19831 : 
iMagorrian fc Tremaind 1 19991 : lAmaro-Seoane et al.l 120041 ') 
where the impact of tidal disruption events is shown to 
be significant and therefore should not be neglected in 
considering the question of what powers the growth of black 
holes. 

There are several mechanism by which stars are 
driven into the loss cone of a black hole. In spherical 
stellar distributions, where the relaxation time T le \ is 
comparable to or smaller t han t he present age of the 
univers^B to (|Freitag et al.l l200ct ). two-body relaxation 
induces a steady change in the angular momentum space 
distribution of stars such that some of them will drift 



1 For the purposes of this study we do not discriminate between 
tp and one Hubble tim e Hq 1 and assume H^ 1 &t = 13.7-10 9 yr 
llKomatsu et al.ll2009l) . 



to very eccentric orbits with pericen tre distances smaller 
than the black hole captu re radius |Frank fc Reesl Il976l : 
iLightman fc Shapiro! Il977t ). In much larger systems like 
the most-massive ellipt ical galaxies which are tho ught to 
be triaxial in shape (|Kormendv fc Bended 119961 s ) . stars 
on box orbits c an cross the central regio n arbitrary close 
to the SMBH (|Binnev fc Tremaind koQgl ) such that they 
become disrupted or s wallowed as a whol e for t he case of a 
very massive SMBH. iMerritt fc Vasilievl (|201Ch concluded 
that the feeding mod e of v ery massive SMBHs, like M87 
|Gebhardt fc Thomas! l2009h . is currently dominated by 
stellar captures. The true rates could be even higher 
since their analysis takes only stellar orbits within the 
black hole influence radius th into account, whereas stars 
within the loss cone but from much further away should 
reach the black hole, too, as long as the critical radius 
?~crit ( a q uantity which is d efined i n Eq. [3t remain s larger 
than r h . iNorman fc Silk! (fl983T): iPoon fc Merrittl (|200ll . 
l2002t) : lMerritt fc Poonl (|2004l ~ iBerczik et al. I (|2006l ) provide 
additional information on the dynamics of SMBHs and 
stellar capture rates in triaxial potentials. 



Observed disruption events ([Ulmerl 19991: Komossa 



2, 
2( 



2002; Halpcrn et al. 2004: Komossa ct al. 2004; Esquei et al 



20081 : iGezari et al.ll2 008: Cap pelluti et a l. 2009; Gczar i et al 



20091 : iKomossa et al.1 120091 : Ivan Velzen et al.1 12010| and ref- 
erences therein), support the view that tidal disruptions 
contribute to the growth history of SMBHs. To which 
magnitude this is the case is a major aspect of this study. 



The paper is organized as follows. In §[5] we will shortly 
explain the concept by which stars are driven by angular mo- 
mentum diffusion into the "loss cone" of the SMBH. This 
formalism is applied to spherical stellar distributions with 
arbitrary slope parameters of the density profile. § [3] de- 
scribes the NBODY6 code that we used. We will specify 
the scale-free models and motivate the very large set of per- 
formed simulations required to infer the disruption rate of 
stars by SMBHs in the nuclei of galaxies. The results will be 
given in § [3] while more detailed information regarding the 
dynamics of the simulations will be part of § [5] In § [6] the 
procedure how to scale the obtained results to realistic as- 
trophysical systems as well as the number of expected tidal 
disruption events will be specified. A critical discussion of 
potential error sources in § [7] is followed by a summary of 
our main findings in § [S] 



2 THEORY 



iFrank fc Reesl (|l976l ) calculate that massive black holes can 
grow not only by swallowing stars which lose their energy 
via dynamical relaxation, but also by swallowing stars on 
very eccentric orbits i.e. stars with low angular momentum. 
The change of the stellar distribution in angular momentum 
space is expected to progress much faster than the change in 
energy space for stars within the critical radius r cr it- Consid- 
ering stars with very low angular momentum and pericentre 
distances smaller than the capture radius of the black hole, 
the velocity vectors of these stars must be aligned very nar- 
rowly. This narrow region is known as the loss cone. It only 
reflects the geometry in velocity space and is characterized 
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by the loss cone angle 6\ c whose symmetry axis is directed 
towards the position of the black hole. For distances below 
the influence radius ru of the black hole where the velocity 
profile follows a Keplarian one (oc r -0,5 ), 9\ c is given by 



i\ c oc 



2/Y 



3r 



(1) 



according to lFrank fe Rees! (|l976l > ) . For r > th a slightly dif- 
ferent expression has to be used. At the moment we leave it 
undefined if the stars are disrupted before entering the hori- 
zon of the black hole or if they are swallowed as a whole. A 
general capture radius r c a p can be specified for all purposes 
|Novikov fc Frolovlll989l ; iBinnev fc Tremaineiliooi ; see also 
Appendix |A)) . In perfectly spherical potentials i.e. poten- 
tials where no torques from anisotropic matter distributions 
can induce an additional supply of stars, all stars on loss 
cone orbits would be consumed within one orbital time scale 
tcross- However, dynamical relaxation between stars causes a 
steady change of the stellar distribution in angular momen- 
tum space and therefore changes in t he velocity vectors b y 
small amounts #Diff per crossing time (| Frank fc Rees|[l976h : 



7Diff OC 



tr 



(2) 



The critical radius r cr i t which is the characteristic distance 
to the black hole where the drift in the velocity vector of a 
star due to dynamical relaxation within one crossing time is 
of the same order as 6\ c is therefore defined as: 



?Diff 



(3) 



Assuming a number density profil^fl n(r) = nor Q within 
r'crit < th and considering only equal mass stars, an expres- 
sion for the critical radius 



■it oc 



Mi 



M 2 n 



(4) 



is obtained by insert ing Eq.[T]and Eq.[2]into Eq.[3l Spitze r's 
relaxation formula (jSpitzer fc HarnJl 19581 : ISpitzerlll987h is 
used for the relaxation time t le \. The Coulomb logarithm is 
neglected. 

The stellar cap t ure ra te can be derived by using eq. 17 
from lFrank fc Rees! (|l97# l: 



(5) 



For a density profile n(r < r cr it) = nor™ the stellar disrup- 
tion rate C is obtained by substituting r = r cr i t : 



C <xG*M. 2 r 



Mi n 



(6) 



For very massive black holes the critical radius becomes 



2 The parameter no can be substituted by no 
more common number density n c at the influence radius rjj- 



We replace v(r) oc 



larger than the influence radiu s of the black hole an d Eq. [4] 
must be modified according to lFrank fc Rees] (|l976h : 



(7) 



r cr it oc (r cap r H no) x + a . 

We assume the velocity dispersion to be a 2 oc 
Gn KLr 2+a and use the same formalism (Eq. [5| to derive 
Eq. [7] The capture rate for r cr i t > th becomes: 



GM(r) 



C oc r C!ip r„ng G* M t 2 (r cap r H n ) w+-> 



(8) 



For simplicity we adopt the number density profile to be 
n(r < r C rit) = nor a and thus assume a to remain constant. 
Real galactic nuclei with SMBHs more massive than 10 7 Mq 
can deviate from pure power law profiles at radii r < r C rit, 
whereas the inner density profiles of large elliptical galaxies 
are nev ertheless well appr oximated by simple power law 
profiles (|Truiillo et al.ll2004 ). Eq.|8]is valid for -3 < a < -1. 

These equations which are derived from the very gen- 
eral a ngular momentum diffusion concept of iFrank fc Reesl 
(| 19761 ), will lose their applicability for systems where the 
stellar phase space is not well-occupied with sufficient 
amounts of low angular momentum stars. Gaps in the phase 
space distribution, for example carved out by binary-SMBH 
evolution must first be repopulated, whereas the relaxation 
driven refilling process may take longer than one Hubble 
time Hq for large elliptical galaxies. Hence for these sys- 
tems the two-body relaxation driven capt ure rate will be 
strongly suppressed (|Merritt fc Wan3l2005l ) . 



3 DESCRIPTION OF THE N-BODY MODELS 

In the following sections the computations and results will 
be specified. We make use of conventional N-body units 
(|Heggie fc Mathieu|[l98r3 ) . For readers being unexperienced 
with these units, a very short overview is given below. 



3.1 N-body units 

The set of N-body units is defined by 



G = M = 1 



(9) 



where G is the gravitational constant and M the total mass. 
If the system is gravitationally bound and in virial equilib- 
rium with r v i r = 1 then the total energy E, which is the 
sum of the kinetic and potential energies of all particles, is 
E = — 1. N-body timescales which are used as the time base 



for the computations are defined to be 



2V2 



Here t c 



tcr. 



is the crossing time of the particles at r = r vlI . The half 
mass i.e. half light radius r e for a constant ^--ratio is usu- 
ally scaled to equal r e — r v ir = 1. It is of the order of kpc 
scales in physical units for real elliptical galaxies. For exam- 
ple one N-body timescale would correspond to t = 7 ■ 10 6 yr 
in physical units for a spherical bulge or galaxy of Sersic type 
(n=4, see § 13. 2\ with a half light radius r e fts 0.65 kpc and 
total stellar mass M = 1O 9 M . In § ED and Appendix [B] 
the detailed procedure how the computational results are 
transformed from N-body to physical units is given. 
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Figure 1. Scale-free density profiles of different Sersic models. 



3.2 Generation of the models 

The observed surface brightness profiles I n (r) of bulges and 
elliptical g alaxies are well approximated by the following 
Sersic law (|Sersidll968l : ICaon et alJll993l ): 
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"/ r N 
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f-1 





Here n is the Sersic index. It represents the strength of 
light concentration towards the center. The parameter 
Te = I(r e ) specifies the surface brightness at the corre- 
sponding half light rad ius r e , whereas b n is a scaling factor 
|Ciotti fc Bertiri Il999h . The 2D density profile can also 
be reconstructed from the measured surface brightness 
profile using an appropriate mass-to-light ratio, which for 
our purposes is assumed to be constant along the radial 
distance to the center of the galaxy. 

In order to study the environmental impacts of massive 
black holes, an unaltered and original Sersic density profile 
is chosen for the initial state of the models. These N-body 
mo dels are set up us ing the same method as described 
in iHilker et al.l (|2007l ). First the 2D Sersic models are 
deprojected into 3D density distributions using Abel's 
integral equation. From the 3D density profile p(r), the 
potential, <f>(r), and mean mass within radius r, M(< r), 
can be deduced. For a non-rotating, spherical system 
with an isotropic velocity distribution, the distribution 
function f(H) is ergodic, where H is the Hamiltonian 
i.e. the total energy of the system. The radial velocity 
distri bution is then deriv e d from f(H) by using eq. 4.46a 
from iBinnev fc Tremainel |2008l ). The actual positions as 
well as velocities of the N-body particles are distributed 
correspondingly in 3D. The program is modified by adding 
a 1/r-potential of the black hole of mass M, = 0.01 in 



N-body units |Heggie fc Mathieul ll986t Fl. This step is 
necessary because otherwise the velocities of particles close 
to the black hole in the N-body computations would be too 
low and the system out of equilibrium. The cut off radius 
for the models is chosen to be 20 times the half light/mass 
radius. 



3.3 NBODY6 numerical dynamics software 

The up to date version of NBODY6 (|Aarsethlll99"9l . l2003h 
with Graphical Processing Unit (GPU) support is used for 
the direct N-body integrations. A black hole is added by a 
SMBH particle of mass M. = 0.01. It is implemented into 
the models at the center of mass while being initially at rest. 
Particles which fall below the limit of the capture radius 
r™ are removed from the simulations while their masses 
are added to the SMBH particle. The capture radius rj!™ 
remains constant- To ensure correct dynamics, the SMBH 
particle receives the center of mass velocity after the capture 
event. 



3.4 The need for a large set of simulations 

In order to extrapolate many scale-free models to astro- 
physical systems which contain some orders of magnitudes 
more stars than are possible to be simulated with direct N- 
body integration methods on modern GPUs, the relaxation 
driven effects in angular momentum and energy space as 
well as every other N dependent systematic effect (see § 15.11 
& § 1 5 . 2 p must be determined. This can be achieved by 
simulating models with different numbers of particles but 
otherwise identical physical parameters. In doing so sev- 
eral particle models following a Sersic n = 4 density pro- 
file are generated. It is desirable to simulate these models 
for as many different black hole configurations as possible 
in order to use the formalism in § 16.11 for the extrapola- 
tion to the black hole of interest, hence increasing the com- 
putational effort considerably. The masses of the particles 
rrii — iV -1 are always scaled to ensure m » = 1 in N-body 
units (|Heggie fc Mathieulll986h . N=15xlk, 15x2k, 10x5k, 
5x10k, 5x25k, 2x50 k, 2x75k and one model containing 
each 100 k, 150 k, 250 k and 500 k particles are generated 
and simulated. All these models are simulated forward in 
time up to ^= crossing times at the virial radius r v i r = 1 
i.e. 100 N-body timescales for three different black hole cap- 
ture radii = 2,4,8 • 10~ 7 . Energy values and relative 

energy errors \AE\ = ~^^p~~T~^^ are evaluated directly 



4 The differences between a 1/r-potential and a realistic 
Schwarzschild or Kerr black hole potential are completely insignif- 
icant for distances of several hundred r cap away from the black 
hole. This is typically the distance where the innermost particles 
are located. 

5 In reality the capture radius would change as well as the to- 
tal number of capture events. However in order to simplify our 
extrapolation formalism to realistic galaxies and due to the fact 
that the mass gain of the black hole within T = 100 NBODY 
timescales is limited to the order of a few percent, it is assumed 
to be constant. 
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Figure 2. The capture rates per one N-body time unit for the three different black hole capture radii, evaluated from the total amount 
of swallowed particles within the timespan of T = 50, 75 & 100 N-body time units. These values are best fitted by the power law function 
C(N) = aN b , here N refers to the total number of simulated particles. 



with the NBODY6 software and controlled every new N- 
body timescale. The relative energy errors usually not ex- 
ceeded values of \AE\ — 1CP 8 — 1CP 4 . A few models had to 
be discarded afterwards as they suffered from repetitive en- 
ergy errors in excess of \AE\ = 10 -2 . To guarantee unbiased 
capture rates we also discarded models in which the posi- 
tion of the black hole was offset by a distance d > 0.1 from 
the density center of the particle distribution. The statistical 
significance of the numerous low N models is increased by 
simulating as many models as possible. The required time 
for the computations of all simulations exceeds a timespan 
of seven months on five modern GPUs. 



4 RESULTS 

In Fig. [2] the number of particles being swallowed by the 
black hole is plotted against the total number of particles. 
This is done for each black hole capture radius r™ = 
2,4,8 ■ 10 -7 . Moreover the total number of captured par- 
ticles within T = 50, 75, 100 N-body integration times is 
divided by these values to obtain the capture rate per N- 
body timescale. The number of captures averaged over all 
runs are then approximated by a power law function 

C(N) = aN b (11) 

with the help of the Marquardt-Levenberg minimization 
method and independently by a grid scanning algorithm 
minimizing the Chi-square error statistics. The free pa- 
rameters a and b have to be positive real numbers while 
the boundary condition C(N)\n=o = requires no off- 
set. To reduce the correlation between the parameters 
a and b to zero, we normalize the power law function 
C(N) = a' (N/N L ) b during fitting. The denominator N L 
refers to the logarithmic mean. The resulting effect can 
be seen in Fig. [3] These uncorrelated valuefl are used for 

6 To simplify the extrapolation formalism, the rcnormalized con- 
stant of proportionality a 1 and its error is afterwards transformed 



the extrapolation to realistic values. The justification for 
using a power law approximation for the capture rate C(N) 
from the simulations comes from Eq. [6] when replacing 
no = Npo and M* — N^ 1 . Poisson square root errors y/ N c 
are assumed for all values and N c is the total number of 
captured particles. The results can be found in Table [1] 
Additionally the reduced Chi-Square values Xn an d the 
X 2 -probability function Q{p,X 2 ) are calculated in order 
to test the validity of a power law approximation for the 
capture rate. Given the values in Table [TJ the hypothesis of 
a power law function seems to be a reasonable assumption. 
However, for the determination of the error values of 
parameters a, b the square root errors \/N c are rescaled 
slightly by the values ^/Xm from Table [T] to obtain %n = 1- 
Otherwise the quoted e rror values would be underestimated 
for the case of Xm > 1 (|Press et al.lll992t fl. 

The advantage of numerical simulations over analytical 
expressions like Eq.|6]are given in the ability to take dynam- 
ical aspects like cusp formation, dynamical heating (§ I5.1|l 
and a wandering SMBHs (§ 15. 2[) into account. These depend 
strongly on time and on the total number of particles 
and may influence the capture rate C(N). The predicted 
power-law index of Eq.[6]is therefore not expected to exactly 
match the value obtained from the simulations. In fact Eq. [6] 

would only predict C(N) oc N^*+°? « N 0M for a w -1.5 
compared to C(N) oc JV ' 83 from the computations. The 
difference is caused by stronger dynamical evolution and 
cluster heating in low number particle simulations accom- 
panied by a decrease in the total number of particles falling 
into the black hole. Models containing many more particles 
have much smoother potentials and relaxation driven effects 
(notably cluster heating) need longer to influence (decrease) 
the capture rate (Fig. [4]). Consequently the exponent b of 
the power law function which approximates the number of 

back to a = This does not affect the correlation coefficient 

N l 

p = between a, b. 
7 Chapter 15.1 
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1.01 
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4 


• io- 


7 


14.46 ± 2.09 
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0.245 




8 
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7 


18.15 ± 2.51 


0.833 ±0.012 


1.41 


0.022 


100 


2 


• io- 


7 


8.73 ±1.33 


0.831 ±0.013 


0.83 


0.788 
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• io- 
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10.98 ± 1.49 


0.841 ±0.012 


1.12 


0.255 




8 


■ io- 


7 


14.41 ± 1.97 


0.845 ±0.012 


1.57 


0.005 



Table 1. Fit parameters of the power law approximation (Eq. lllH 
for the simulated Sersic n = 4 models. The black hole capture 
radii and timescales T are given in N-body units. Xn = X 2 /M 
corresponds to the reduced Chi-Square values, /i are the degrees 
of freedom and Q = r(0.5/x, 0.5x 2 ) the x 2_ P r °t , ability function 
which estimates the likelihood of the power law fit. 



captured particles of the total set of simulations becomes 
larger than expected from Eq. [15] These dynamical processes 
are reflected by the values of a, b at different timescales. The 
constant of proportionality a decreases in time, whereas 
the slope parameter 6 is consistent with a small increase 
from b « 0.80 at time T = 50 up to b « 0.83 at time 
T — 75. Thus the exponent of the power law function which 
approximates the capture rates becomes slightly larger, 
whereas the constant of proportionality decreases. Moreover 
the T = 50 values may still be influenced by initial condi- 
tions. There are minor changes in b from T — 75 to T = 100. 

For the purpose of this study the rate C is assumed, 
within the statistical uncertainty, to remain unchanged when 
extrapolated to larger values of N. This assumption can only 
hold if the phase space is already well occupied with suffi- 
cient amounts of low angular momentum stars. This is a 
necessary condition for the steady diffusion process of stars 
into the loss cone. The capture rates of the Sersic n — 4 
models are found to be maximal at the beginning of the 
simulations in contrast to Sersic n = 2 models with their 
much shallower density profiles (Fig. [4]). This demonstrates 
the above assumption to be credible, at least for galactic 
nuclei containing SMBHs less massive than 10 7 Mq. In such 
galaxies the diffusively refill of any small gap with radius 
''gap << th would a nyway occur on a timescale shorter 
than a Hubble time (|Merrittl 120051 ). The observed strong 
JV dependence (b=0.83) may become irrelevant or absent 
for black holes more massive than 1O 7 M0, especially if they 
have core profiles. For these systems the loss cone refilling 
timescale T re fin ~ df c T le i, becomes very long. Once the ini- 
tially filled loss cone becomes emptied within a few crossing 
times, the capture rate C would stagnate at insignificant 
values as long as there is no re-population mechanism more 
efficient than angul ar momentum diffusion (|Merritt fc Wangl 
l2005l ; lMerrittll2005h . 

This effect can be illustrated by simulating Sersic 
n = 2 models. These have a slower dynamical evolution, a 
different cusp and cluster heating timescale and a reduced 
population of low angular momentum stars compared to 
the Sersic n = 4 models. In this way, qualitative limitations 



0.9 

0.85 
0.8 
0.9 

0.85 
0.8 
0.9 

0.85 
0.8 

0.75 



";. r cap= 2e - 7 ' NreP 1 

\T= 100 

1 o — 

2 a * 

3 a 


r cap =2e-7, N ref =N L 


.../cap= 4e - 7 > N re pl 

; X,, " 


ft 

r cap= 4e - 7 ' N ref=N L 


r cap= 8e - 7 ' Nrer 1 


r cap =8e-7, N ref =N L .... 

ft 



5e-5 le-4 1.5e-4 2e-4 1 1.2 1.4 1.6 1. 

a 



2 2.2 



Figure 3. The error ellipses for the models after T = 100 before 
(left) and after (right) renormalization. The shape of the error 
ellipses becomes nearly circular which proofs the parameters to 
be uncorrclatcd. 
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Figure 4. Time evolution of the capture rates for Sersic n = 4 
& n = 2 models. The statistical significance of the latter ones is 
increased by averaging over three simulations. 



on the number of capture events for core-type galaxies 
with shallow central density profiles (Fig. [TJ can be ob- 
tained. Even though the extended outer profiles of the 
most-massive elliptical galaxies are conform with a large 
Sersic index n, the 'depleted' core-type central regions (this 
is where the relevant black hole physics take place) are 
more similar in their appearance to the shallow centers 
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of low n model^f]. A strongly reduced disruption rate in 
comparison to the Sersic n = 4 models is evident in these 
computations. The enlarged radius of influence r« and 
therefore the difference in the extrapolation formalism to 
realistic galaxies can not compensate these differences. 
Moreover in the largest simulated Sersic n — 2 models, the 
capture rate stagnate first around insignificant values. It 
starts increasing (Fig. [4j afterwards, accompanied by the 
relaxation driven formation of a cusp and a population 
of stars with sufficiently low angular momentum. If we 
assume this behaviour to persist unchanged up to even 
larger numbers of particles, i.e. to large core-type galaxies 
where no cusps can form on timescales shorter than Hq 1 , 
these numerical findin gs confirm analytical predictions 
l|Wang &: Merrittl I2004T ) in a qualitative way. The capture 
rate of stars in large core-type galaxies is very low, as long 
as the diffusive refill of the angular momentum space with a 
sufficient number of stars, i.e. the cusp formation timescale, 
takes longer than a Hubble time. 

See also §[?]& Appendix [Cl for more details on this topic. 

Finally the here performed simulations of the Sersic 
n = 4 models strongly support the scenario of lFrank fc Reesj 
(| 19761 ) in which stars are driven into SMBHs via diffusion 
in angular momentum space and not only by diffusion in 
energy space. From the most elementary considerations of 
energy diffusion and by assuming the two-body relaxation 
time to be T rc i oc lr ^ N ^ , one would expect C(N) = ^jj- oc 
TV • T~, oc ln(7V). Such a small increase of the capture rate 
with N is incompatible with our results. 



5 DYNAMICS & SCALING ISSUES 

The capture rate is influenced by several dynamical pro- 
cesses which are described below. 



5.1 Cusp formation and cluster expansion 

The process of relaxation st rongly influences the d ynamics 
of stars around a SMBH. In lBahcall fc Wold (| 19761 ) the re- 
laxation driven evolution of the stellar density profile near 
a SMBH is determined. It is found that the energy which 
some stars loose through near encounters is balanced by an 
outgoing flux of energy if the slope of the density profile 
is a = —1.75. The required time to form such an equilib- 
rium density B&W profile strongly depends on the relax- 
ation time which becomes larg er the smoother a gravita- 
tional potential is (|Spitzerlll987r ). 



The a = —1.75 profile is compatible with the present 
N-body models only up to TV = 25 — 50 k where relaxation is 



8 This is also one reason which complicates the discrimination be- 
tween 'true' cores, formed by the dynamical evolution of massive 
binary black holes and those in which only the outer envelopes 
are modified by near encounters. In the latter case the outer pro- 
file ext rapolated to inwards radii w ould suggest the existence of 
a core l|Hopkins fc HernquisfclOf ). The binary black hole mech- 
anism may also be acc ompanied by other processes lowering the 
centr al stellar density l|Merritt fc Vasiliev|[2uTol : ISchawinski et alj 
12006ft . 
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Figure 5. Time evolution of the central slope parameter a plot- 
ted for the 50 k, 150 k and 250 k models. The linear trend (solid 
line) is only drawn for the 150 k and 250 k models, while the 
first one (fixed black hole) was simulated forward in time up to 
T = 200. 



strongest and the statistical scatter is large. The TV > 50 k 
models, which allow a more precise measurement of a, are 
found to be in the developing stage towards more cuspy pro- 
files. In Fig. [5] the time dependent central slope parameter 
a within r = 0.004 is plotted for some models. The radius 
r is chosen to be 20% smaller than the time and model- 
averaged black hole influence radiu^f] in order to ensure that 
the slope parameter is not determined for radii larger than 
th at the beginning of the simulation when the mass and 



9 See Appcndix[B]for information regarding the determination of 
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Figure 6. Mass profiles of two models. The thin dashed black line represents the gradient of the B&W profile while the thick dashed 
black line (only drawn for T=0) displays the unaltered Sersic n = 4 model. The first error on a corresponds to the fitting error while 
the second one to the statistical error inferred from Monte Carlo simulations. The profiles and thus a are evaluated for radii r < 0.004. 
For more informations sec the text below. 



influence radius of the black hole is smallest. In order to 
obtain the central slope parameter a, it is inappropriate to 
calculate the density profile p(r) oc r a from given shells of 
thickness Ar and densities p(r + Ar). Unfilled shells, espe- 
cially in low N models, would strongly bias the determina- 
tion. In order to circumvent this difficulty, the cumulative 
mass function M(r) oc r 13 oc J Q r r 2 por a dr is calculated and 
the density slope parameter a = ft — 3 (equating coeffi- 



cients) is determined from the measured ft. This approach 
is tested by Monte-Carlo simulations in which several thou- 
sand models of particles following a p oc r 1-1 ' 76 distribution 
are realized. For each of these models the central slope pa- 
rameter a within r — 0.004 is calculated. The models are 
scaled such that the number of particles within r = 0.004 
is equal (within the statistical scatter) to those of the 25 k, 
50 k, 75 k, 150 k and 250 k simulations. The standard devi- 
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Figure 7. A comparison between the time evolution of several 
Lagrange radii for three different simulated models. As expected 
from theory, the Lagrange radii evolve faster to larger values in 
simulations containing fewer particles. The fluctuations are sta- 
tistical in nature. The position of the black hole is used as the 
reference center. 

ation a from the obtained normal distributiorPl of central 
slope parameters is then taken as a reasonable estimate for 
the statistical error in addition to the one obtained from the 
fit itself. In Fig. [6] the time evolution of the mass profiles of 
two models are plotted. 

In order to estimate the dependence of a wandering 
black hole on cusp formation processes and finally the 
capture rate, simulations of fixed black holes are desirable. 
Such simulations are realized by making use of a modified 
NBODY1 code (see £15.31 for more details regarding the 
capture rates). Within the large statistical errors, no 
significant difference in the density profiles between the free 
floating and fixed black hole is identified for the 50 k model. 
This is not an unexpected finding since the most bound 
particles, which are also the particles with the highest 
probability of being captured, are expected to follow the 
motion of the black hole. However a rigorous statistical 
evaluation is beyond the scope of this study. 

While the capture rate is increased by cusp formation, 
dynamical heating counteracts by reducing the central 
density. The cluster starts to expand by decreasing the 
absolute value of its binding energy due to increasingly 
more strongly bound particles which are losing energy by 



10 Actually very small particle numbers within r = 0.004 bias the 
power-law density-approximation and the distribution of central 
slope parameters becomes asymmetric with a tail towards very 
large values. This may partially account for some extreme out- 
liers especially in low N models, whereas for larger models the 
distribution becomes more symmetric and the expectation values 
fi center around a = —1.75. 



relaxation. These particles, which may finally be swallowed 
by the black hole, are transferring their kinetic energy 
to other particles. This heating is illustrated by the time 
evolution of the Lagrange radii (Fig. 0. As a consequence 
the capture rate is expected to depend strongly on the 
density profile close to the black hole (Eq. [SJ. 

In reality mass segregation of heavier bodies being 
relev ant for multi-mass s yst ems (Alexander fc HopmarJ 



20091; iBaumgardt et all 



— jpr_- .. 

2004bT iMorrisI Il993l ; 

Preto fc Amaro-Seoan stellar collisions 

1 Bailey fc Daviesl Il999l ; iDaleet al l 120091). a s ignificant 
fraction of primordial binary stars jHopmanll2009l ). torques 
from aniso tropic matter dis tributions acting as massive 
erturbers JPerets et al. 20071 ) . star formation by gas inflow 

and the possible presence of 
120061 ) would complicate the 



201C 



Hopkins fc Quataertl 
IMBHs jBaumgardtet all 
dynamics of stars close to a SMBH even more. These effects 
are also expected to accelerate the dynamical evolution and 
to enhance the number of stellar disruption events. Newly 
formed stars may replace those lost by tidal disruptions 
while tidal torques from IMBHs or a second SMBH are 
expected to refill the loss cone efficiently. Recoiled black 
holes should also e nforce a burst of stellar disruptions 
(|Stone fc Loebll2010h . In nature the relaxation driven B&W 
cusp formation takes very long and is expected to exceed 
one Hu bble time H^ 1 for black hole masses larger than 
1O 7 M (|Freitag et alJl200Sh . 



5.2 Wandering black hole 

In the simulations the SMBH particle responds to the 
interaction with other particles which causes the SMBH to 
wander. This might affect the fo rmation of a density cus p 
and influence th e capt ure rate (|Baumgardt et ail l2004al ) . 
IChatteriee et all (|2002i ) gives a very detailed overview of 
the relevant forces acting on a SMBH. They are summarized 
below. 

The here performed simulations differ only in two ways 
from the N-body simulations done by IChatteriee et al] 
(|2002l ). The black hole is allowed to swallow particles and 
the forces are unsoftened. The SMBH moves around the 
common center of mass due to the gravitational interaction 
with particles bound to it, whereas unbound particles are 
forcing the black hole to wander in a way which resembles 
the Brownian motion of molecules. The latter process is the 
dominant contribution to the wandering of the black hole 
(see Fig. lU . 

The situation is now complicated by the possible 
occurrence of violent three body encounters, e.g the inter- 
action between the black hole, a strongly bound particle 
in orbit around it and another one. Recoil events force the 
black hole and its surrounding particles to move outwards. 
The mass fraction, is usually orders of magnitudes 

larger in any performed N-body simulation than it is in 
a realistic nucleus of a galaxy. And, because the recoil 

and for 



effect becomes stronger for a larger fraction 
smaller capture radii r^™, wandering of the SMBH in the 



M. 
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Figure 8. The 100 binned (Ad = 0.001) x, y, z-positions of the 
SMBH particle with a capture radius r" a ™ = 4 ■ 10 -7 for the 
500 k, 250 k and 75 k model. In the last row the sum of these 
values is plotted and approximated by a normal distribution. The 
probability distributions are well approximated by a Gaussian 
underlining the character of the Brownian motion. The relevant 
length scale d is given in units of the virial radius r vl[ = r e = 1. 
The SMBH particle in the 75 k model experienced a minor kick 
during the integrations. 



Figure 9. In the upper two figures the binned specific energy — 

-7 



distributions of the 50 k, r" a ™ 



: 4- 10- 7 and 500k, r=™ = 4-10" 



models are plotted. The lower diagrams depict the ratio of cap- 
tured particles to total number of particles within the given en- 
ergy bins. Evidently only the particles with the most negative 
energy i.e., the most strongly bound particles are accreted as ex- 
pected from theory. The upper and lower (black) lines represent 
the error uncertainties. 



N-body models is expected to modify the processes leading 
to the formation of a cusp. Consequently, the simulated 
capture rate may become affected. If the recoil kick of 
the SMBH particle is strong enough to eject it out of the 
density center or even from the whole cluster, the capture 
rate would drop significantly. This is expected, due to 
obvious reasons, to happen more likely in simulations with 
low particle numbers. As a consequence the extrapolated 
N-dependent capture rate would be strongly biased and the 
best fitted slope parameter, b, may be too large. Therefore 
the actual position of the SMBH particle is compared to the 
density center of the matter distribution for every simulated 
model and at every new N-body time unit. The black hole 
particle is not considered in the calculation of the density 
center which is determ ined by the method described in 
ICasertano fc HutJ (|l985l ). If the position of the black hole 
and the density center are offset from each other by d — 0.1 
in N-body units, the simulation is removed and replaced 
by a different one. In nearly all simulations this offset is 
smaller than 10 -3 — 1CP 2 . This guarantees that the results 
are not biased by displaced black holes in the low N models. 

But even by removing those few models where "un- 
natural" kicks and displaced black holes are observed, the 
wandering of the black hole itself might affect the capture 
rate. The wandering radius can be determined by the 



standard deviation of the normal distribution (Fig. [8). It 
is found to be comparable in size to the influence radius 
th = 0.005 (for the 250 k model) and becomes gradually 
smaller for larger particle numbers i.e. smaller mass frac- 
tions ^. 

A first clue about the degree to which the wandering 
black hole affects the results can be obtained by a closer 
look at the energies of accreted particles. If only particles 
are swallowed which are strongly bound i.e. have the most 
negative energies, the effect of Brownian motion on the 
capture rate is expected to be rather small, since the cloud 
of strongly bound particles moves together with the black 
hole. In Fig.|9]the initial energy distribution for two models 
is shown. Also plotted is the fraction of the accreted parti- 
cles to the total number of particles within a given energy 
bin. Evidently only the most strongly bound particles are 
captured. If the energy E = -^P- + Q.5mv^ + 0.5M.vl u 
of the particle of mass m and black hole is negative, shortly 
before it enters the capture radius and is removed, the 
particle is gravitationally bound to the SMBH. In our 
models the vast majority of particles are gravitationally 
bound to the black hole, e.g. the fraction of bound particles 
centers around 100% in the low-N models and 85 - 95% in 
the largest-N models. 

We therefore conclude that a wandering black hole 
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does not bias the capture rate in a way that would make it 
unrealistic when extrapolated to real IMBHs and SMBHs. 
The performed simulations automatically contain the 
gradual change in the number of accreted particles which 
are influenced by the wandering of the black hole. Our 
largest N-computations already approach realistic IMBHs 
embedded in globular clusters. To resolve all doubts that 
the steep dependence on N of the capture rate, C oc TV 83 , 
is not caused by the systematics of the wandering black 
hole, especially in low N models, direct N-body simulations 
with fixed black holes (§ 15. 3[) are performed. For complete- 
ness it should also be mentioned that the N-body models 
include two additional effects: (i) A restoring force which 
arises between the black hole and the overall potential of 
the stellar distribution, especially if it has a cuspy density 
center, and (ii) a dynamical frictional force when the black 
hole p asse s through the cloud o f particles (| Chandrasekharl 
Il943allbllc3 ; IChatteriee et~alll2002h . 
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5.3 Fixed black hole 

Simulations with a fixed black holes are realized by us- 
ing NBODYI. Unfortunately it is impossible in NBODY6 
to fix the SMBH particle to a specific location while si- 
multaneously using all of its computational benefits. On 
the other hand the usage of an independent N-body soft- 
ware implementation reduces the possibility of systematic 
errors. The NBODYI simulati ons are performed on special- 
purpose, GRAPE-6A boards (|Fukushige et alj|2005l ) at the 
stellar Populations and Dynamics Research Group in Bonn. 
The black hole is mimicked by an (unsoftened) external - 
potential which is directly implemented into the code. Parti- 
cles which cross the capture radius are removed, while their 
masses are added to the mass of the black hole. To circum- 
vent collisions between field particles, a small softening pa- 
rameter e = 10~ is used. Additionally some strongly bound 
particles around the external potential are erased artificially 
(the number corresponds to roughly 30% of the total num- 
ber of "true" capture events) in order to prevent gradual 
slow downs, large energy errors and/or the complete crash 
of the simulations. The energies of all removed particles are 
handled carefully to ensure a correct energy output. Due to 
these limitations and the much smaller sample of simulated 
models, the NBODYI computations are not used for the 
extrapolation to realistic galaxies but only for a rough com- 
parison to the much more advanced NBODY6 simulations. 
In Fig. [10] the results are plotted. Despite the large simplifi- 
cations of the NBODYI computations, the power law index 
b of the capture rate, C oc N b , agrees, within the statistical 
uncertainties, closely with the index obtained with the much 
more sophisticated NBODY6 simulations with free moving 
SMBHs. As a consequence a (strongly) wandering black hole 
particle does not bias the low N results in a way which would 
be dangerous when extrapolating these to astrophysical sys- 
tems harboring many more stars than particles in our simu- 
lations. Of course this behaviour may change for initial black 
hole masses different from the one M. (f = 0) = 0.01 used in 
these computations. 



Figure 10. Results of the simulations with a fixed black hole. 
Strongly bound particles which needed to be removed artificially 
to prevent slow downs or a computational crash are not considered 
in the evaluation of the slope. Note the excellent agreement with 
the results obtained for the more realistic NBODY6 computations 
(Fig. 0. 



6 DISCUSSION 

6.1 Scaling to realistic galaxies 

The so-far presented results must be scaled to astrophys- 
ical systems in order to infer the rates at which stars are 
disrupted by central, supermassive black holes. From the 
following relation 



Th Isim 



Th I astro 



(12) 



which must be necessarily fulfilled, the capture radii, r- cap , 
for the corresponding black holes of interest must be deter- 
mined. In order to sca le to astrophysical systems, we use the 
AI. - a relation from lSchulze fc Gebhardtl (|201lh . 



M. 



1.51 



, 200km s" 1 

and the expression for the radius of influence, 

GM. 

r H = — 5-, 

G 

to calculate rg for a SMBH of given mass, 
13.1 



Ms J [P4 



(13) 



(14) 



(15) 



Here Mg corresponds to 10 s Mq and for reasons of com- 
putational feasibility we neglected the intrinsic scatter of 
the Mm — <J relation. This is useful when dealing with av- 
eraged quantities like the impact of stellar disruptions for 
the growth history of the majority of SMBHs. Some studies 
may instead be interested in individual systems and the ex- 
trapolation formalism can easily be replaced by direct mea- 
surements of M.,th and a instead of using the values from 
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the Mm — a relation. This also holds for the choice of the 
relevant tidal disruption radius, 

reap = gn 3 , (16) 

where g is a parameter depending on the stellar polytrope, 
mass and spin of the black hole as well as the trajectory 
of the star. Black holes below 10 7 Mf) and solar like 
stars are well approximated by g ~ 1 jKochanekl 1 19921 ) . A 
more detailed discussion of r cap can be found in AppendixlAl 

The relevant astrophysical timescale is ob- 
tained through the computation of the crossing time 
tcr(rff) = alrlj ) ^ ^he influence radius m of the black 
hole in comparison with that of our numerical integra- 
tions. The number of disruption events within the given 
timescale is obtained from the derived capture rate 
C(N,rf™) = a(rf™)N b . Here N refers to the total number 
of (real) stars with the averaged mass M* in the bulge 
component or whole elliptical galaxy. It is assumed to be 
N = 10 j/' f * in accordance with our numerical integrations, 
whereas a(rf^) is extrapolated to the black hole mass of 
interest by using the T — 100 values for the parameters <£3- 
The parameter b is assumed to be unrelated to the capture 
radius and is hence taken to be constant at b = 0.83 ± 0.01. 
In fact Eq. [S] predicts the slope parameter b to be unrelated 
to the capture radius. Nevertheless a minor change in b 
towards smaller vales of r^p 1 cannot be rejected given the 
b values of Table \T\ at T — 100. This might be explained 
by a combination of timing issues, simplified assumptions 
of our analytical approach or is purely statistical in nature. 
Therefore the parameter 6 is extrapolated (by linear and 
power law regressions) down to the required values of r™ 
in order to test its impact on the capture rates. The impact 
is found to be moderate because r™ nas t° be extrapolated 
down to rflp = (0.06 - 0.07) • 10~ 7 (depending on the used 
M. — cr-relation) for the largest black hole with 10 7 Mq. 
While the capture rate would be unaffected for the least- 
massive black holes, it would drop by a factor of 2 for the 
most-massive ones. Increasing uncertainties of these values 
due to the propagation of error analysis strongly overlaps 
with those of fixed b. For the purposes of this study we 
therefore assume the parameter b to be independent of r™ 
and refer the reader to § [7] for a more critical discussion on 
that topic as well as of the improvements left for future work. 

Finally for an individual galactic nucleus hosting a 
SMBH of mass M,, with a radius of influence th, velocity 
dispersion o(r — m), capture radius r cap and a stellar pop- 
ulation with the mean mass M*, a very general expression 
for the capture rate inferred from the numerical integration 



11 The capture rates C(N) from the numerical computations are 

normalized to one N-body time unit i.e — y= crossing time at the 

2 v 2 

virial radius r vlr = 1 and must be scaled down to one crossing 
time at the influence radius of the black hole in order to become 
synchronized with the astrophysical timescale t CTOSS . 

12 At least three different capture radii must be simulated to al- 
low for non linear extrapolation of the parameter a(r^™)- This is 
required for the extrapolation to different black hole sizes/masses. 



can be obtained by applying Eq. 1141 

s 0.951 

r< n oooci / JU * 1 -1.363 0.363 0.363 

Castro = 0.00061 I — I r H r* g a (17) 

The validity of Eg. 1 171 covers the parameter range of IMBHs 
as well as SMBHs up to M. ~ 10 7 M© . In the following 
section we explicitly make use of the M. — a relation and 
assume only solar like stars as well as g — 1. 

6.2 Disruption rates of IMBHs & SMBHs 

By applying the extrapolation formalism from section £16.11 
the integrations yield the following expression for the cap- 
ture rate of real astrophysical galaxies: 

/ M \ °' 446 

C(M.) = 6.29.10- 8 (jLy [yr-]. (18) 

For comparison the results are also extrapolated accord- 
ing to an older versi on of the M, — a relation from 
iFerrarese fc Fordl (|2005l ) to illustrate the dependence of the 
capture rate on systematic black hole mass determinations: 

/ M \ 353 

C(M.) = 3.54 • lO" 7 (^J [ yr -i]. (19) 

These results holds for nonrotating isotropic galaxies 
or globular clusters with cuspy inner density profiles 
p(r) oc r a , where the density power law index is a ~ —1.5. 
Eq. [18] and [19] should not be applied to black holes with 
masses larger than 1O 7 M0. The uncertainties correspond 
to about 50% of the values (see Table [2}. The interested 
reader is referred to Appendix [B] for a much more detailed 
description of how the numerical results are extrapolated to 
realistic galaxies. The astrophysical disruption rates of stars 
(including the statistical uncertainties) for some exemplary 
black holes are summarized in Table [2] We also calculate 
the disruption rates for I MBHs in order to compa re them 
with previous simulations (|Baumgardt et ah I l2004al ). 

The expected number of tidal disruption events in galac- 
tic nuclei containing black holes of 10 6 to 10 7 Mq inferred 
from the numerical integrati ons are in good agreeme nt with 
recent optical based surveys ijvan Velzen et al.ll2010h . While 
their study yields the rate for tidal flares per galaxy to be 
C = 3(11,) • 10 _5 yr _1 , the results obtained from the present 
simulations give C = 3.0(±1.4) - 8.3(±4.2) • lO'V" 1 for 
black holes in the mass range 10 6 to 1O 7 M . 

The simulations also offer some clues about the growth 
of IMBHs and SMBHs in the lower mass range. We observe 
only a modest impact of the black hole mass on the capture 
rate. For a mass range over four orders of magnitude, the 
capture rate increases only by a factor 25-60 depending 
on the used scaling relation. The relaxation driven growth 
of massive black holes by stellar disruptions is thus only 
important for IMBHs and SMBHs up to several 10 5 M Q . 
IMBHs should easily double their mass within a few Gyr 
in per fect agreement with earlier studies (jBaumgardt et aD 
l2004al ) . Much more massive black holes must have grown 
by different processes rather than the relaxation driven 
infall of stars, i n goo d agreement with the findings of 
lYu fc TremaTnel (|2002T ) and gas accretion and feedback 
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Schulze fc Gebhardt (2011) 



M.(1O 6 M ) 


c(io-V _1 ) 




L(ergs 


- 1 ) 


U.UUl 


n i a. A- n f\« 
U.l*± m u.ud 


n 11 -L n 

U.ll I u.uo 


O Q _l_ 1 "7 


1U 


U.Ul 


n -i-ni7 


u.oy zc u. io 


i i _i_ n ^ 
1 . ± m u.o 


1U 


0.05 


0.8 ±0.4 


0.9 ±0.4 


2.2 ± 1.0 


10 40 


0.1 


1.1 ±0.5 


1.4 ±0.7 


3.0 ± 1.4 


10 40 


0.25 


1.6 ±0.8 


2.3 ± 1.1 


4.5 ± 2.2 


10 u 


0.5 


2.2 ± 1.0 


3.4 ± 1.6 


6.2 ±3.0 


10 40 


1 


3.0 ± 1.4 


4.9 ±2.4 


8.4 ±4.1 


1Q 40 


2 


4.0 ±2.0 


7.3 ±3.6 


1.1 ±0.6 


10 41 


4 


5.5 ±2.7 


11 ± 5 


1.6 ±0.8 


10 41 


10 


8.3 ±4.2 


18 ±9 


2.4 ± 1.2 


10 41 




Ferrarese & Ford C2005) 






A/.(1O 6 M ) 


c(io-V _1 ) 




L(ergs 


- 1 ) 


U.UUl 


0.40 ±0.17 


0.04 ±0.02 


l.i m u.o 


in 4 " 

1U 


U.Ul 


0.90 ± 0.40 


0.16 ±0.07 




1U 


0.05 


1.6 ±0.7 


0.46 ±0.21 


4.5 ±2.1 


10 40 


0.1 


2.0 ±0.9 


0.72 ±0.34 


5.8 ±2.7 


10 40 


0.25 


2.8 ± 1.3 


1.3 ±0.6 


8.0 ±3.8 


1Q 40 


0.5 


3.6 ± 1.7 


2.0 ± 1.0 


1.0 ±0.5 


10 41 


1 


4.6 ±2.2 


3.2 ± 1.5 


1.3 ±0.6 


10 41 


2 


5.8 ±2.8 


5.0 ±2.5 


1.7 ±0.8 


10 41 


1 


7.4 ±3.7 


7.9 ±3.9 


2.1 ± 1.0 


10 41 


10 


10.3 ±5.2 


14 ±7 


2.9 ± 1.5 


10 41 



Table 2. The expected number of stellar disruption events C for 
solar like stars by supermassive black holes up to Mm < 10 7 Mq. 
For comparison our numerical results are extrapolated according 
to an older version of the M. — cr-relation jF crrarcsc fc Fordl2005f l 
and the most recent one llSchulze fc Gebhardtll201ll) . Within a 
factor of two they agree with each other. T^d is the time needed 
to double the initial mass of the black hole in units of the Hubble 
time Hq 1 . Only one half of the stellar mass is assumed to become 
accreted by the black hole (Rccs 1988). Finally the time averaged 
mean luminosity L = O.beCM-QC 2 of these black holes is calcu- 
lated by assuming the efficiency parameter of matter to energy 
conversion to be e = 0.1. The motivation behind is to compare 
these energies with potentially detectable left overs of relativis- 
tic outflows which may become dep osited into the surrounding 
medium after tidal disrupt i on events |C rocker fc Aharoniarl201ll ; 
iGiannios fc Metzgerlfeplll ; | van Velzen et alTl201lT) . However the 
deposited energy strongly depends on the formation rate of rel- 
ati vistic jet out flows and may be significantly overestimated by 
us (Bower 2011). Nevertheless these deposited energies might be 
relevant for studies aiming to make a robust detection of dark 
matter annihilation signals in galactic bulges, dwarf galaxies or 
globular clusters hosting a central black hole. These results have 
relevance for galaxies with cuspy density profiles with slope pa- 
rameters a —1.5 within the inner most few pc. 



mode ls (|Silk fc Rees! Il998l ; iFabiarj 1 19991 ; iMurrav et all 

Our findings exclude any relevance for establishing 
the M, — g relation from stellar disruptions in density 
profiles similar to those of the simulations. This is due to 
the relatively small capture rate and hence large doubling 
times (T2D > Hq 1 ) for black holes more massive than 



3 The growth of the very early population of SMBHs may also be 
domi nated by stellar disruptions in isothermal cusps llZhao et al.l 
120021) . See the information in the text below. 



10 Mq. If for example the initial mass of a SMBH is 
strongly under-massive with respect to the M. — a relation, 
the feeding from tidal disruptions events alone might not 
be sufficient enough to bring it close to the observed 
relation for galaxies at z ~ 0. On the other hand if stellar 
disruptions dominate the growth of the least-massive black 
holes there is no obvious reason why these black holes 
should follow the M, — a relation. By now assuming the 
M. — a relation to be established for a primordial gas 
rich globular cluster (or galactic nucleus), which nowadays 
remains in isolation and without gas to drive new star 
formation, the resulting IMBH (or SMBH, at least if it is 
not too massive) should nowadays be more massive than 
expected from the M, — a relation due to subsequent tidal 
disruption events. It is very te mpting to connec t these 
results to the case of oj-Centauri (|Novola et al.|[201ol ). Tidal 
disruption events might therefore have implications for 
the search and existence of IMBHs in globular clusters. 
Of course in order to proof its relevance for IMBHs, the 
use of the M. — a relation in the extrapolation formalism 
from numerical simulations to galactic nuclei (Appendix [B} 
must be replaced by more direct observational data because 
the extrapolated values strongly depend on the validity of 
this scaling relation. Tidal disruption events complicate 
the understanding of the relevant processes which drive 
the evolution of galaxies and their central black holes. 
Especially as the impact of disruption events for the mass 
growth of black holes strongly depends on their initial mass. 

In spite of this it mi ght be interesting to re late these 
findings to a recent study iKormendv et alj (|201ll ) in which 
observational evidence for secular growth processes of black 
holes in disks and pseudobulges is found. The capture 
rate for rotation-supported models like rotating bulges or 
pseudo bulges should be enhanced compared to nonrotating 
models. These o bjects are expected to form f rom rotating 
bar instabilities jKormendv fc Kennicuttll2004l ) and the rel- 
ative velocities between two or more particles are generally 
lower. Therefore two-body relaxation processes would be 
even stronger. 

However the overall picture of black hole growth across 
cosmic times by tidal disruptions might be complicated 
even more due to the dynamical evolution of the density 
profile and a variable fraction of the initial stellar mass 
which finally becomes accreted by the black hole. Our 
conclusions regarding the growth history of IMBHs and 
SMBHs events should only hold for density profiles re- 
sembling those of our simulations and by assuming that 
a fraction of one half (or more) of th e initial ste llar mass 
becomes accreted by the black hole (Rocs 1988|). In fact 
some effects can considerably reduce this fraction and 
complicate the efforts to estimate the significance of tidal 
disruption events for the overall growth history of black 
holes. Recent hydrodynamical simulations suggest that for 
loss cone stars on nearly parabolic orbits, most of the stellar 
matter is ejected within the first orbit and then later on 
due to powerful shocks which may be energetic enough to 
ignite thermo nuclear reactions unbinding large amounts of 
stellar mass (IBrassart fc LuminetJ l200a iGuillochon et al.l 
120091 ). Secondly and especially relevant for black holes in 
the lower mass range, accretion luminosities far in excess of 
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the Eddington limit jStrubbe fc Quataertll201ll ) may blow 
away most of the remaining gas. In the end the growth of 
these black holes due to tidal disruption events may be 
insignificant even for very large capture rates of several 
events per 10 6 yr. 



7 CRITICAL DISCUSSION AND OUTLOOK 
FOR FUTURE WORK 

To the best of our knowledge this study reports for the first 
time the expected tidal disruption rate of stars by SMBHs 
up to 1O 7 M0 obtained by direct N-body integrations. N- 
body computations offer a large amount of advantages over 
analytical studies. They can handle several physical effects 
simultaneously while most analytical studies are forced to 
simplify at least some of the dynamics. On the other hand 
direct N-body integrations aiming to infer astrophysically 
relevant numbers of stellar disruption events are confronted 
by their own limitations and difficulties. In this section we 
will critically review limitations of our own simulations as 
well as improvements and ideas left for future work. 

(a) In Table [2] we calculate among other values the re- 
quired timescale Tid for doubling the mass of a black hole 
of given initial mass. This timescales is computed from the 
total number of captures averaged over 100 N-body time 
units (see Table [TJ. We recommend the reader to regard the 
doubling time T^d only as some reference guide. When ex- 
pressed in physical time, our simulations last only a fraction 
of one Hq 1 (between several 10 7 and one 10 9 years) and 
may not represent much longer time episodes. Moreover we 
assumed one half of the disrupted star to be accreted by 
the black hole. There exist two effects that can reduce the 
amount of stellar matter which finally becomes swallowed 
by the black hole. First, if the tidal stripping occurs from 
a nearly parabolic orbit, hydrodynamical simulations sug- 
gest one half of its m ass to be lost within its first path 
l|Guillochon et al1l2009l ') and large quantities of the remain- 
ing mass to be bl own away by shocks and thermonuclear 
reactions later on jBrassart fc Lum inct 200§1|). Second, very 
small black holes might temporarily generate luminosities 
far in excess of the Eddington limit (|Strubbe fc Quataertl 
l201lf ) and most of the remaining matter may finally be 
blown away instead of being swallowed by the black hole. 
This would invalidate our conclusions regarding the growth 
history of small black holes where we assumed one half of the 
stellar mass to be accreted. Nevertheless the inferred cap- 
ture rate should be valid for all galaxies or stellar clusters 
with density profiles comparable to our simulated ones. 

(b) With current generations of CPUs it is unthinkable 
to simulate galaxy models with realistic numbers of stars 
with direct N-body integration methods. The only way to 
obtain stellar disruption rates for SMBHs in the centers 
of galaxies is to simulate as many models as possible to 
infer all relevant N-dependent systematics affecting this 
rate. Afterwards the results can be extrapolated. However 
it is important not to do this for only one given black 
hole capture radius but for many black hole configurations. 
Therefore all these simulations must be repeated for several 
capture radii r^S in order to extrapolate them according 
to the formalism in § 16.11 to the black hole of interest. 



We calculate the capture rate for black hole masses in the 
range 1O 3_7 M0. Due to the highly nonlinear Eq. 1131 we 
had to extrapolate parameter r™ from Table [T] down to 
r™ « 0.07 ■ 10~ 7 . The usage of three different black hole 
capture radii is thus the minimal requirement to obtain 
useful values under the assumption that the parameter 
a ( r cap) from Eg. II II follows a power law distributional with 
positive parameters and no offset. 

There is no question that future studies must redo these 
simulations for different capture radii to constrain a(r^p) 
even more precisely. However this is a very time consuming 
task. The complete set of our Sersic n — 4 simulations took 
more than seven months to compute on five modern GPUs. 
Despite the large amount of needed computing power, 
direct integration methods like NBODY6 may exceed their 
limitations when the capture radius falls significantly below 
10 -7 in N-body units, especially if the mass of the black 
hole particle is of the order of one percent or more of the 
total masf^l. In addition to that the statistics may worsen 
(due to a limited number of capture events) and must be 
balanced by even more simulations. 

In our computations no severe r^ 1 ™ dependence of the pa- 
rameter b is evident, in accordance with theoretical consid- 
erations (Eq.[6lfc[8lF|- Therefore we assumed it to be con- 
stant. However we cannot exclude per se any deviation at 
very small capture radii. A systematic decrease in the pa- 
rameter b for even smaller values of r^™ would only reduce 
the tidal disruption events of the more massive SMBHs in 
our sample. We plan to tackle this problem as well as to 
constrain the parameters a(r^L™) and b even more precisely 
in the future. 

(c) In this study effects from General Relativity are ne- 
glected. The relevant tidal disruption radius of a SMBH for 
solar like stars is several times larger than its Schwarzschild 
radius and relativistic effects should become strongly 
suppressed for radii r » r s . This makes our assumption 
of neglecting GR credible. Nevertheless a fully relativistic 
treatment of a black hole potential yields a deeper gravita- 
tional potential than a purely Newtonian one, thus being 
more attractive for compact bodies like stars to be captured 
by the SMBH. On the other hand particle scattering by 
a relativistic potential may result in stronger deflection, 
perhaps powerful enough to reject some stars from the 
immediate vicinity of the black hole thereby decreasing the 
capture rate. The next generation of N-body integrators 
is expec ted to be sop histicated enough to address these 
aspects (|Aarsethll2007l ). 

Rotating black hole spacetimes should be considered, too. 
The cross section of a realistic black hole strongly depends 
on its mass and spin parameter j = jf-, where J. is the 
angular momentum of the black hole. The likelihood for a 
particle to become swallowed by the black hole depends on 



14 a ( r cLp) is specified in Eq.lB3l 

15 Private communication with Sverre Aarseth. 

16 According to the theory of angular momentum diffusion, pa- 
rameter b only depend on the slope parameter of the density pro- 
file. 
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the spin parameter, its trajectory and angular momentum. 
Particles are more likely captured if they counter rotate the 
black hole because in this direction the effective capture ra- 
dius is enlarged. It is not unreasonable to conclude that a 
rotating black hole embedded inside a nonrotating spherical 
distribution of stars will lose some of its angular momen- 
tum. On the one hand j decreases when M. becomes larger, 
on the other hand counter rotating particles are more likely 
captured. This would lower J. and thus the spin parameter 
j. If tidal disruption events really contribute a significant 
amount of mass to a specific population of black holes, it 
should also affect their spin values in a way which might 
deviate from the predictions of gas accretion models. 

(d) A crucial quantity for extrapolating our numerical re- 
sults to astrophysical systems is the black hole radius of in- 
fluence th- For its evaluation we use the kinematic determi- 
nation (Appendix [Bj. We observe this radius to be roughly 
five to six times smaller than the dynamical radius r g . This 
is the radius at which the mass in stars/particles equals the 
mass of the black hole. If interested readers plan to rescale 
our models by replacing the M, — a relation by directly mea- 
sured data of th for some galaxies, it is very important that 
they also use the same influence radii as the ones used in 
our simulations and not the dynamical radii. 

(e) The capture rate from our numerical results should 
not be applied to SMBHs above 10 7 M Q . The refill of the loss 
cone takes a timespan of the order T rc fin « #f c T re i. The refill 
of the loss cone is much faster in N-body integrations than 
in nature, since the potentials are more cuspy and relaxation 
times are shorter than in reality. Therefore our simulations 
have only relevance for galactic nuclei where the loss cone 
refilling times are much shorter than Hq 1 . In Appendix [C] 
we show that Trefiil << Hq for a black hole with a mass 
of 1O 7 M . This becomes also evident from Eq. [6]&[8] For 
black holes significantly more massive than 1O 7 M0, i.e with 
very large particle numbers and very smooth potentials, the 
critical radius becomes much larger than the influence ra- 
dius of the black hole and Eq.[6]has to be replaced by Eq.[8] 
The latter one predicts a different behavior for C(N) such 
that the numerically found capture rate should not be ex- 
trapolated to black holes in excess of 10 Mq. By inserting 
the relevant values from or computational findings to the 
systems of interest, Eq. [5] predicts the critical radius not 
to exceed the influence radius for black holes less massive 
than 1O 7 M0, thus showing our simulations to be governed 
by processes r cr i t < rjj. 

(f) One could even criticize the black hole mass M,(t = 
0) = 0.01 used for our numerical computations to be too 
high as the black hole m ass fraction in realistic galaxies is 
a factor of a few smaller (|Magorrian et alii 19981 ). Neverthe- 
less most of the relevant dynamics happens at distances of 
the order of the influence radius th whereas we use the 
radius of influence for the extrapolation to realistic galax- 
ies. The choice of M.(t = 0) = 0.01 is therefore not ex- 
pected to change the capture rate significantly. In this con- 
text the usage of different capture radii instead of differ- 
ent initial masses M,(t = 0) for the extrapolation to the 
wide set of astrophysical SMBHs should be justified, too. 
The strict relation between mass and capture radius of a 
black hole (Eq. IA1[) enables variation of the latter one while 
keeping the former one constant in the scale-free N-body 
simulations. The great advantage of this strategy is given 



in equal black hole influence radii, crossing times, cusp for- 
mation timescales etc. simplifying the extrapolation formal- 
ism considerably. The same holds true for the overall Sersic 
n = 4 profiles. Not every outer bulge component or ellip- 
tical galaxy profile resembles that of a Sersic n = 4 i.e. de 
Vaucouleurs profile. Mostly relevant for the direct number 
of capture events is the density profile close to th ■ For relax- 
ation times smaller than one Hq 1 the formation of a cusp 
(up to a = —1.75) is expected. Such a gradual change of the 
density profile is also found in the numerical simulations. 
Hence our simulations cover a large space of isotropic, non- 
rotating density profiles for black hole masses up to 10 7 Mq . 

(g) We only treat single-mass systems while galactic cores 
are known to be multiple-mass systems featuring additional 
processes like mass segregation, star formation, binary evo- 
lution, torques from anisotropic matter distributions, reso- 
nances etc. Stellar remnants like neutron stars would not be 
disrupted outside the event horizon and could probe much 
deeper potentials than solar like stars, thus complicating 
the gravitational dynamics and making relativistic correc- 
tion terms inescapable. They would also disappear without 
any visible counterpart when finally captured. 

(h) Finally our numerical simulations should only be 
regarded as a first (very) limited approach to a systematical 
scan of capture rates in galaxies. It would be important 
to extent these studies by simulating the same models for 
even smaller capture radii rtext and longer timescales in 
order to reduce the need of extrapolation. It would be 
important to take into account rotating and triaxial stellar 
density profiles around the SMBH and to decrease the still 
rather large uncertainties. Direct N-body simulations of 
isothermal p(r) = 2 J Gr ^ spheres, which might represent 
the initial phases of elliptical galaxies an d bulg es best, 
should be performed as well. IZhao et al.l (|2002t ) found 
evidence for strong black hole growth in isothermal cusps. 
By assuming r- cap oc r s = 2 a" * , rewriting Eq. [5] to 
C(r) cc p{r)r 2 o6h by using r„ = ^ and 9f c = 

for very massive SMBHs, one obtains C(r) oc ■ (J-f-) 2 ■ 
Under the assumption that the capture rate is dominated 
by stars from th, the r dependence cancels out and the 
final mass of the black hole is M m (tf) = J+Zq* C{r = 

r H )At » 10 8 M Q • (—S^) 8 ^.). This relation is 

indeed i n very close agr eement to the observed M, — a 
relation ()Zhao et al.ll20oj ). Therefore stellar captures might 
contribute significantly to the growth of SMBHs in the 
past, especially if the loss cone refill is enhanced by mergers 
and/or triaxial stellar distributions. 



Despite some of the details stated above, the here re- 
ported simulations represent (the first) systematic estimate 
for the capture rate by SMBHs of stars in galaxies with 
cuspy inner density profiles. This work should be followed up 
by simulating different capture radii r^£p as well as density 
profiles, taking relativistic correction terms into account, 
and by trying to find ways to infer the numbers of disrup- 
tion/capture events for SMBHs with mass > 10 7 Mq. 
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8 CONCLUSION 

We performed direct N-body simulations to obtain the num- 
ber of disruption events of stars by SMBHs which are pre- 
sumed to exist in the centers of most galaxies. A modified 
NBODY6 code was used. All computations were processed 
by several CPUs over several months integration time. The 
initial density profiles of the models were chosen to follow 
nonrotating isotropic Sersic n = 4 profiles. We calculated 
numerous models with different particle numbers but oth- 
erwise equal physical parameters in order to ensure good 
statistics. This is required because all systematic effects de- 
pending on the total number of particles must be specified in 
order to extrapolate the simulations to realistic galaxies by 
using the formalism presented in § 16.11 The rates at which 
stars are captured are found to be nearly independent of 
the mass of the black hole. Thus only the growth over cos- 
mic times of IMBHs and of the least massive SMBHs may 
be dominated by stellar disruptions. The expected tidal dis- 
ruption rate is a few events every 10 5 years per galaxy for 
black holes in the mass range up to 10 7 Mg. The feeding by 
stars from density profiles similar to the ones computed here 
bears no implications for establishing scaling relations be- 
tween very massive black holes and their host galaxies. This 
is in agreement with conventional gas accretion/feedback 
models. On the other hand the growth history of the least 
massive black holes might be governed by more than one 
feeding mode (gas and star accretion). This might have im- 
plications for the search and existence of potential IMBHs in 
globular clusters and minor galaxies. Assuming these scal- 
ing relations (e.g. the M, — a relation) to be established 
shortly after their primordial gas rich phase billions of years 
ago, the nuclear black holes would continue their growth by 
the subsequent disruption of stars. Depending on the ini- 
tial conditions, the black hole masses could nowadays lie 
well above the predicted values of the M. — a relation as 
long as the globular cluster remains in isolation^ 7 ! On the 
other hand the continuous monitoring and search for tidal 
disruption events in globular clusters (e.g. in the Virgo Clus- 
ter) should constrain the fraction of those clusters hosting a 
central IMBH. By assuming 25000 globular clusters with a 
central black hole in the mass range M. = 10 3 - 10 4 M Q m 
the Virgo Cluster of galaxies, there should be one disruption 
event every 10 — 25 years. Finally the performed computa- 
tions indicate that the growth history of IMBHs and low 
mass SMBHs is diverse and not only governed by one pro- 
cess, i.e gas accretion. However it needs to be pointed out 
that there exist effects which might reduce the fraction of 
stellar matter which finally becomes accreted by the black 
hole. We a ssumed one half of a captured star's mass to be 
swallowed (|Reeslll988h . whereas a smaller fraction would re- 
sult in even slower growth rates. Thus our conclusions re- 
garding the growth history may change if small black holes 
gather only tiny fractions of the total initial stellar mass. 
Future studies can use the reported capture rate C(M.) to 
deduce more realistic growth rates M(M.) by taking more 
appropriate values for the fraction of accreted matter into 
account. It would also be interesting to extend these studies 



The relevant velocity dispersion a should therefore not in- 
crease. 



to the most-massive black holes as well as constraining the 
capture rate for different profiles. 
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the required capture radius 



for a black hole of mass 



M m must be obtained by using astronomical observations of 
indi vidual galaxies or by making use of the M. — a relation 
from lSchulze fe Gebhardtl (|201lh . If in the near future much 
larger samples of measured SMBH masses allow for more ac- 
curate values, it will be no problem to implement them into 
this formalism. By combining Eq. IB 1 1 with the disruption 
radius r cap = gr^^^f-^ and the expression for the radius 

0.54 

[pc] which is derived from 



of influence th ~ 13.1 



(mI) 



the M. 



a scaling relation, r^p 1 follows: 



APPENDIX A: THE TIDAL 
DISRUPTION/CAPTURE RADIUS 

The disruption radius at which a star is torn apart by tidal 
forces such that roughly one half of its matter will become 
accreted by the black hole, is a function of the mass and 
spin of the black hole as well as the trajectory, internal 
structure, size and mass of the star. A star is disrupted 
outside the event hori zon if the black hole mass is smaller 
than a certain limit (|Lai et al.l Il994l ; iBinnev fc Tremaine! 
120081 ) . For typical solar- type stars with masses Af* « 1M© 
and radii r* « 1R© the mass of the black hole must be 
smaller than M. < 10 s Mq to disrupt the star before reach- 
ing the event horizon. A stro ngly spinning black hole dra - 
matically alters the situation (|lvanov fc Chernvakovall2006f ). 
Sufficiently massive black holes swallow stars as a whole. The 
General Theory of Relativity predicts the radius where a star 
is doomed to enter a very massive black hole to be larger 
than the actual Schwarzschild-radius r a (|Novikov fc Frolovl 
Il989h . Stars on initial Keplarian orbits coming from infinity 
with pericentre distances q < 4r s will be swallowed by the 
black hole as long as (— ) <S 1. The particles in the N- 
body simulations do not come from infinity but their speed 
at the apocentre distance is much lower than the correspond- 
ing speed of light and the capture radius of r cap = 4r a seems 
to be the most natural and best approximation for the be- 
haviour of a realistic (extremely massive > 1O 8 M0) black 
hole. This approximation also holds for the bound particles 
around the black hole which are most likely swallowed. The 
ratio ( " ap ° ) 2 for apocentre distances of 10 -4 — 10~ 2 is al- 
ways much smaller than one. For our purposes the capture 
radius can finally be defined as: 



gr, 

BGM. 

A 



(mI) 



M. < 10 8 M G 
M. > 10 8 M G 



(Al) 



The parameter g, which is of the order of one, dep e nds o n 
many paramete rs and can be taken from iKochanekl jl992tl : 
lLai et all l|l994f ): llvanov fc Chernvakoval (|2006t ). 



APPENDIX B: EXTRAPOLATION 

In the following part we give a more detailed description 
of the formalism by which the here obtained capture rates 
(Table [TJ can be scaled up to realistic bulges of galaxies or 
elliptical galaxies. 



(a) From the relation 



Th I sim 



T H I astro 



(Bl) 



4g ■ 10" 



M, 

Ms 



(B2) 



It specifies the required capture radius in the scale-free N- 
body integrations for the astrophysical black hole of interest. 
Afterwards the function a(rc' a ™) must be evaluated from the 
values in Table [T] 



a(rip) = 0.023(±0.006) (r£j}) 



0.363(±0.020) 



(B3) 



yields a reasonable approximation 18 ! for the extrapolation of 
the parameter a from Eq. QT]to anv desired rf™. For the 
purposes of this paper the slope parameter b = 0.83 is as- 
sumed to be independent of r^L™. As already mentioned in 
§ 16.11 the parameter g accounts for the stellar model and 
mass of the black hole. It is of the order of one jKochanekl 
1 19921 : lLai et al.|[l994l ; llvanov fc Chernvakoval200ri ) . For sim- 
plicity we use g — 1 which is a reasonable assumption for 
nonrotating black holes less massive than M, — 10 7 Mq and 
solar mass stars. Eq. IB2I assumes all stars to be disrupted 
before entering the horizon. 

(b) The dynamical timescale t a i m of the N-body particles 
inside the sphere of influence rg has to be calculated ac- 
cording to tsim = ,' 2th n ~ 0.008. It is used as a reference 
for timing issues when compared to the relevant astrophys- 
ical timescales t. To ease the extrapolation of the numeri- 
cal results to astrophysical systems, we compute the time 
averaged influence radius th- Representative for all mod- 
els we calculate th and t s i m from the 25 k, 50 k, 75 k, 150 k 
and 250 k models. For the calculation of the radius of in- 
fluence we bin the particles in cylindrical shells of thick- 
ness Ar = 0.001 and measure for each configuration the 
one dimensional velocity dispersion (line of sight velocity) 



in order to obtain cr(r)g im . Here Ni is the num- 
ber of particles within each configuration. We choose the 
line of sight axis to be parallel to the z-axis. Afterwards 

°bh,i = """f]v~^ ' 7") i s calculated for each cylindrical 

shell to obtain cr(r) 2 h , here n = \J x\ + y'f + z'f. The factor 
3 in the denominator is used for the normalization to the 
relevant line of sight velocity inside the isotropic distribu- 
tion. The radius of influence rg is then calculated to be the 
radius at which , .°j m = 2. We note that in N-body units 
G = 1. The position of the black hole is used as the refer- 
ence center and the mass gain of the black hole is taken into 
account. For the time averaged influence radius and velocity 



18 Q = 0.89 without rescaling Xn = 1- Afterwards the uncertain- 
ties are taken directly from the covariancc matrix. Renormaliza- 
tion induces the errors to be uncorrelated to each other. 
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dispersion we obtain ~ 0.005 and a(r = th) ~ 1.26. The 
black hole influence radius is 5 — 6 times smaller than the 
dynamical radius. 

(c) Subsequently the astrophysical dynamical timescale 
tcr(rtf) = 2l£Lj of the matter distribution within the 

CI v " ' ct I astro 

influence radius of the astrophysical galaxy must be com- 
puted for the black hole of given mass by using th ~ 



13.1 (jfc) [pc] a nd a 



200 



[kms 



from 



ISchulze fe Gebhardtl (|201lh . 

(d) The number of stars N in the astrophysical galaxy 
must be specified. For simplicity we assume all stars to have 
the same mass (M+) = IMq. A coarse estimate for the num- 
ber of stars can be computed by: 

100A/. 



N = 



(B4) 



The choice of (M*) = IMq depends on the stellar mass func- 
tion and seems to be a reasonable assumptio n for galactic 
nucle i where mass segregation is importa nt (|Freitag et alj 
120061 ; lKroupall200ll ; iLSckmann et al.ll2010h . The factor 100 
accounts for the fraction of bulge mass to black hole mass 
in accordance with our simulations. 

(e) Finally the disruption rate of stars by massive black 
holes can be evaluated. In a first step the numerically in- 
ferred number of captures C(N) per N-body time unit (Ta- 
ble [TJ must be normalized to the relevant crossing time 
isim = 0.008 (in N-body time units) at the influence radius of 
the black hole. This dimensionless number must afterwards 
be synchronized with the relevant timescale t cr (rjf) of the 
astrophysical galaxy. Consequently C(N) ■ 4im has to be di- 
vided by tcr(fir) in order to obtain the number of disrupted 
stars within the desired physical time unit (e.g yr, Myr) for 
the black hole of interest: 



C a; 



0.008 ■ a(r?£)N b 
icr(rff) 



(B5) 



Our extrapolation formalism strongly depends on the M. — a 
relation. More accurate and numerous black hole measure- 
ments will improve this relation in the future. Moreover we 
only treat errors from our simulations and neglected the in- 
trinsic scatter of the M. — a relation for simplicity. 




-► r 

Figure CI. Sketch of a typical loss cone problem 



<*r cap <- 



be r c 



M. = 10 7 M o 



, the loss cone angle 8 can be evalu- 



ated from Eq.[T] For the constant of proportio nality / we use 
/ = 2 in accordance with lFrank fc Reesl(| 19761 ). By assuming 
th, Mj, ~ lMfl , the relaxation time to be T ra \ ~ HZ 1 



|Freitag et al.l 120081 ) and estimati ng all other relevant pa - 
rameters from the M. — a relation l|Ferrarese fc Fordll2005l ). 
one obtains the desired result T re fiii ~ 4- lO -6 -!^ 1 << ^o" 1 - 
Even though our assumptions are idealistic and not every 
star fulfills its plunge into the black hole from the critical 
radius r or it) it underlines the extrapolation from our numer- 
ical simulations to realistic cores of galaxies with central 
black holes up to 10 7 Mq to be credible. 



APPENDIX C: LOSS CONE PROBLEMS 

Direct N-body integrations are limited by a maximal com- 
putable number of particles which is orders of magnitudes 
lower compared to particle numbers in the nuclei of as- 
trophysical galaxies. The extrapolation to such astrophys- 
ical settings is thus only possible if the relevant physics do 
not change in between. Loss cone problems (tidal captur- 
ing and/or shrinking binary black holes) require special care 
|Gualandris fc Merrittll2007h . Here we will show that the in- 
equality 

Irefiii = 2 T Ic i << H 1 (CI) 

is fulfilled up to SMBHs of order 10 7 Mq and hence our re- 
sult, C oc A^ ' 83 , should yield realistic values when extrapo- 
lated to such black holes. For black holes much more mas- 
sive, the situation might change. By assuming the radius 
r at which particles can enter loss cone orbits without be- 
ing scattered away through interactions with other stars to 
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